#进行题目的分类，即每个分类的题目考察同类的属性 只与Q矩阵有关
#number是初始分类个数的限制，题目类别对应Q矩阵的行数，自定义函数主要得出各个分类每个类别下的具体的题目
##进化终极版本
match_rows_end <- function(response, number) {
  J <- ncol(response) # 题目个数
  colnames(response) <- paste0("Q", 1:J)
  
  same_matrix <- matrix(NA, J, J)
  rownames(same_matrix) <- colnames(same_matrix) <- paste0("Q", 1:J)
  
  # 计算相似度矩阵（相同作答数）
  for (i in 1:J) {
    for (j in 1:J) {
      same_matrix[i, j] <- sum(response[, i] == response[, j])
    }
  }
  diag(same_matrix) <- 0
  
  # 找出每一列的最相似题目
  max_positions <- sapply(1:J, function(i) which.max(same_matrix[, i]))
  df <- data.frame(
    item = paste0("Q", 1:J),
    max_item = sprintf("Q%d", max_positions)
  )
  
  matched_groups <- list()
  processed <- logical(nrow(df))
  
  # Step 1: 初步聚类（基于相似连接）
  for (current_row_index in 1:nrow(df)) {
    if (processed[current_row_index]) next
    current_group <- as.character(df[current_row_index, ])
    processed[current_row_index] <- TRUE
    
    repeat {
      new_overlap_found <- FALSE
      for (i in 1:nrow(df)) {
        if (!processed[i] && any(current_group %in% as.character(df[i, ]))) {
          current_group <- unique(c(current_group, as.character(df[i, ])))
          processed[i] <- TRUE
          new_overlap_found <- TRUE
        }
      }
      if (!new_overlap_found) break
    }
    matched_groups <- append(matched_groups, list(current_group))
  }
  
  # Step 2: 若初步分类数 > 要求数，执行合并
  while (length(matched_groups) > number) {
    
    # 计算组中心
    group_centers <- lapply(matched_groups, function(group) {
      rowMeans(response[, group, drop = FALSE])
    })
    
    n_groups <- length(matched_groups)
    distance_ma <- matrix(NA, nrow = n_groups, ncol = n_groups)
    rownames(distance_ma) <- colnames(distance_ma) <- paste0("G", 1:n_groups)
    
    # 计算欧式距离矩阵
    for (i in 1:(n_groups - 1)) {
      for (j in (i + 1):n_groups) {
        distance_ma[i, j] <- sqrt(sum((group_centers[[i]] - group_centers[[j]])^2))
      }
    }
    
    # 找最小距离
    min_index <- which(distance_ma == min(distance_ma, na.rm = TRUE), arr.ind = TRUE)[1, ]
    i_min <- min_index[1]
    j_min <- min_index[2]
    
    # 合并两组
    new_group <- unique(c(matched_groups[[i_min]], matched_groups[[j_min]]))
    
    # 更新分组列表
    matched_groups <- matched_groups[-c(i_min, j_min)]
    matched_groups <- append(matched_groups, list(new_group))
  }
  
  return(matched_groups)
}
# 模拟作答数据DINA------------------------------------------------------------------------------------------------
# 输入被试的属性掌握模式 Q矩阵 s g 失误猜测参数
DINA_model <- function(attribute_profile,Q_matrix,item_s,item_g) {
 # browser()
  N <- nrow(attribute_profile)
  J <- nrow(Q_matrix)
  # 被试属性掌握模式 统一为矩阵
  attribute_profile <- as.matrix(attribute_profile)
  Q_matrix <- as.matrix(Q_matrix)
  # num_people <- nrow(attribute_profile)
  # num_items <- nrow(Q_matrix)
  probabilities <- matrix(NA, nrow = N, ncol = J)
  for (i in 1:N) {
    for (j in 1:J) {
      required_attributes <- which(Q_matrix[j,] == 1)
      # ideal_response <- attribute_profile[i, required_attributes] == 1
      #ideal_response <- (attribute_profile[i,])[which(Q_matrix[j,] == 1)] == 1
      ideal_response_vector <- attribute_profile[i, required_attributes]
      if (all(ideal_response_vector == 1)) {
        probability <- 1 - (item_s[j])
      }else {
        probability <- item_g[j]
      }
      probabilities[i, j] <- probability
    }
  }
  Response <- matrix(NA, nrow = N, ncol = J)
  for (i in 1:N) {
    for (j in 1:J) {
      a <- runif(1)
      Response[i, j] <- ifelse(probabilities[i, j] > a, yes = 1, no = 0)
    }
  }
  return(Response)
}


#输入邻接矩阵，输出可达和属性掌握模式----------------------------------------------------------------------------------------------------
#Ad_matrices为邻接矩阵
ad_Reachability <-  function(Ad_matrices) {
  k <- ncol(Ad_matrices)
  A1 <- Ad_matrices + diag(k)
  B <- (A1 + diag(k) > 0) * 1
  #邻接矩阵加单位矩阵等于初始的可达矩阵
  # 可达矩阵扩张算法 循环结束后B为可达矩阵 
  repeat {
    R <- A1 %*% B  # 矩阵乘法
    R[R > 1] <- 1
    if (identical(R, A1)) {  # 如果R不再变化，则停止迭代
      break
    } else {
      A1 <- R  # 更新A1
    }
  }
  reachability_matrix <- R
  # 扩张算法根据可达矩阵算Qs
  #browser()设置交叉进行布尔值判别
  indicate <- ncol(R)
  for (i in 1:ncol(R)) {
    for (j in 2:ncol(R)) {
      if (i < j) {
        new_col <- as.matrix(as.numeric(R[,i] | R[,j]))
          if (!any(apply(R,2,function(x) all(x == new_col)))) {
            R <-  cbind(R,new_col)
            indicate <- indicate + 1
          }
      }
#有新列产生进行判别，新列加入进行属性掌握模式的更新
    }
  }
  
  # 加入00000 代表什么也没掌握的people
  zero_col <- matrix(data = 0,ncol = 1,nrow = ncol(B))
  zi_matrix <- t(R)
  attribute <- t(cbind(R,zero_col))
  # attribute <- as.matrix(attribute)
  # colnames(attribute) <- paste0("A",1:nrow(attribute))
  return(list(reachability_matrix =reachability_matrix ,test_item = zi_matrix, attribute_matrix = attribute))
}


# 给定单个子矩阵 和 题目分类  形成多种符合条件的Q矩阵----------------------------------

single_Q_matrix <- function(zi_matrix,matched_item,number_items) {
  #browser()
  all_Q_matrices <- list()
  zi_matrix <- as.matrix(zi_matrix)
  number_Q_row <- nrow(zi_matrix)
  number_contribute <- ncol(zi_matrix)
  # 生成所有可能的排序
  assignment_per <- permutations(number_Q_row, number_Q_row)
  # 倒置的Q_matrix 会有l1的阶层种情况
  for (j in 1:nrow(assignment_per)) {
    current_assignment_order <- assignment_per[j, ]
    Q_matrix_contrary <- matrix(data = NA,nrow = number_contribute,ncol = number_items)

    # 对于每一种分配方案
    for (m in 1:length(matched_item)) {
      cols_to_fill <- as.numeric(sub("Q(\\d+)", "\\1", matched_item[[m]])) # 提取题号
      # 取当前分配的属性模式（如第 current_assignment_order[m] 行）
      label_to_use <- zi_matrix[current_assignment_order[m], ]
      Q_matrix_contrary[, cols_to_fill] <- label_to_use
    }
    Q_matrix <- t(Q_matrix_contrary)
    all_Q_matrices[[j]] <- Q_matrix
  }
  return(all_Q_matrices)
}

#  从邻接矩阵的角度来生成符合条件的属性层级
#  先判断是否为有向无环图
is_DAG <- function(adj) {
  #browser()
  # 获取最大节点
  k <- nrow(adj)
  # 统计每个节点的入度，即每列中1的个数，即为该节点的边数
  indeg <- colSums(adj)
  # 找到入度为0的节点 作为拓扑排序的起始点
  queue <- which(indeg == 0)
  # 初始化计数
  visited <- 0
  # 循环处理，只要队列中还有节点就继续循环
  while (length(queue) > 0) {
    # 取出队列中第一个节点，并将其在队列中移除
    node <- queue[1]; queue <- queue[-1]
    visited <- visited + 1
    # 查找子节点 找到所有由当前node直接指向的子节点
    children <- which(adj[node, ] == 1)
    for (j in children) {
      #节点node被成功排序移除 指向j入度减1
      indeg[j] <- indeg[j] - 1
      # 如果子节点入度降为0，说明所有前驱节点排除j变成新起点
      if (indeg[j] == 0) queue <- c(queue, j)
    }
  }
  # 输出结果true 或 false
  visited == k
}
generate_linjie_matrix <- function(k) {
  if (k < 2) stop("k must be >= 2")
  
  # 上三角部分的元素数量（不含对角线）
  n_edges <- k * (k - 1) / 2
  
  unique_structures <- list()
  seen_signatures <- character()
  skipped_all1 <- 0L
  
  total <- 2^n_edges
  
  # 穷举所有上三角可能性
  for (pattern in 0:(total - 1)) {
    # 二进制转向量（长度 n_edges）
    bits <- as.integer(intToBits(pattern))[1:n_edges]
    
    # 初始化邻接矩阵
    adj <- matrix(0, k, k)
    adj[upper.tri(adj)] <- bits  # 只填上三角
    
    # 跳过“全上三角为1”的情况
    upper_sum <- sum(adj[upper.tri(adj)])
    if (upper_sum == k*(k-1)/2) {
      skipped_all1 <- skipped_all1 + 1L
      next
    }
    
    # 由于只考虑上三角，本身就是DAG，不必再检测 is_DAG()
    
    # 可达矩阵签名（判等价结构）
    R <- ad_Reachability(adj)$reachability_matrix
    sig <- paste(R, collapse = "")
    if (sig %in% seen_signatures) next
    
    # 存储唯一结构
    seen_signatures <- c(seen_signatures, sig)
    unique_structures[[length(unique_structures) + 1]] <- adj
  }
  
  return(unique_structures)
}

# 穷举Q矩阵 根据子矩阵即典型项目考核模式 和可达矩阵
genarate_Q_matrix <- function(zi_matrix,matched_item,keda_matrix,number_items) {
  #browser()
  library(gtools)
  library(dplyr)
  all_Q_matrices <- list()
  zi_matrix <- as.matrix(zi_matrix)
  # 典型项目考核的维度 肯定比可达维度多
  number_Q_row <- nrow(zi_matrix)
  # 题目分类的维度 必须抽几种的题目
  number_factor <- length(matched_item)
  # 可达的维度
  number_keda <- nrow(keda_matrix)
  # 如果考核维度大于分类维度 必须考虑可达矩阵 生成的Q矩阵必须包含可达
  # 如果题目维度等于可达维度 则直接从可达中抽
  # 等于的时候 子矩阵就是可达
  # 当小于的话直接从层级层面就筛除了 在此不考虑
  if (number_Q_row >= number_factor && number_factor == number_keda) {
    assignment_per <- permutations(number_keda,number_keda)
    for (j in 1:nrow(assignment_per)) {
      current_assignment_order <- assignment_per[j, ]
      # 行列颠倒 便于存储矩阵
      Q_matrix_contrary <- matrix(data = NA,nrow = number_keda ,ncol = number_items)
      
      # 对于每一种分配方案
      for (m in 1:length(matched_item)) {
        cols_to_fill <- as.numeric(sub("Q(\\d+)", "\\1", matched_item[[m]])) # 提取题号
        # 取当前分配的属性模式（如第 current_assignment_order[m] 行）
        kd_ma <- t(keda_matrix)
        label_to_use <- kd_ma[current_assignment_order[m], ]
        Q_matrix_contrary[, cols_to_fill] <- label_to_use
      }
      Q_matrix <- t(Q_matrix_contrary)
      all_Q_matrices[[j]] <- Q_matrix
    }
    return(all_Q_matrices)
  }else if (number_Q_row >= number_factor && number_factor > number_keda) {
    
    # 只有在该条件下会出现原始矩阵题目改变的情况
    # 如果变化情况有6种 每一种根据顺序不同有24种情况 一共会有6×24种Q矩阵
    result_list <- list()
    # 其实可以不用分情况
    #先确定出来可达矩阵外，在哪个范围内抽取
    zi_df <- as.data.frame(zi_matrix)
    kd_ma <- t(keda_matrix)
    kd_ma <- as.data.frame(kd_ma)
    
    other_item <- anti_join(zi_df,kd_ma,by = paste0("V",1:number_keda))
    # 抽取数量 除可达行的缺失题数
    dis_number <- number_factor - number_keda
    other_id = 1:nrow(other_item )
    # m = 2 代表每次抽取 2 个行号
    row_combinations <- combn(other_id, m =  dis_number)
    extract_rows <- function(col) {
      # 使用行索引提取对应的行数据
      return(other_item[col, ])
    }
    
    all_row_combinations_list <- apply(X = row_combinations,2 , FUN = extract_rows)
    # 组成新的题目矩阵   再进行Q顺序排除与上述逻辑一样
    new_item <- list()
    for (i in 1:length(all_row_combinations_list)) {
      new_item[[i]] <- rbind(kd_ma,all_row_combinations_list[[i]])
    }
    # 对于每个组合都有在不同位置排列的形成Q矩阵
    assignment_per <- permutations(number_factor,number_factor)
    for (i in 1:length(new_item)) {
      for (j in 1:nrow(assignment_per)) {
        current_assignment_order <- assignment_per[j, ]
        # 行列颠倒 便于存储矩阵
        Q_matrix_contrary <- matrix(data = NA,nrow = number_keda ,ncol = number_items)
        for (m in 1:length(matched_item)) {
          cols_to_fill <- as.numeric(sub("Q(\\d+)", "\\1", matched_item[[m]])) # 提取题号
          # 取当前分配的属性模式（如第 current_assignment_order[m] 行）
          now_new_item <- new_item[[i]]
          label_to_use <- now_new_item[current_assignment_order[m], ]
          if (is.list(label_to_use)) {
            if (length(label_to_use) == 0) {
              warning("sampled_items 列表为空，无法生成 Q 矩阵行数。")
              return(label_to_use) # 返回空结果
            }
            label_to_use <- do.call(rbind, label_to_use)
          }
          Q_matrix_contrary[, cols_to_fill] <- label_to_use
        }
        Q_matrix <- t(Q_matrix_contrary)
        all_Q_matrices[[j]] <- Q_matrix
      }
      result_list[[i]] <- all_Q_matrices
    }
  }
  return( result_list)
}

# 海明距离
HDD <- function(x, IMP, Qmatrix) {
 # browser()
  N <- nrow(x)
  NI <- nrow(IMP)
  J <- nrow(Qmatrix)
  
  if (!is.matrix(x)) x <- as.matrix(x)
  if (!is.matrix(IMP)) IMP <- as.matrix(IMP)
  
    
  # 标准 IRP 计算逻辑 (DINA/DCM 模型): IRP = 1 if all required attributes are mastered
  check_mastery <- IMP %*% t(Qmatrix) 
  IRP <- 1 * (check_mastery == matrix(
    rowSums(Qmatrix),
    nrow = NI,
    ncol = J,
    byrow = TRUE
    ))
    if (!is.matrix(IRP)) IRP <- as.matrix(IRP)

  
  # --- 3. 计算海明距离并分配 ID ---
  
  # distances 维度：N (参与者数) x NI (IMP数)
  distances <- matrix(NA, N, NI)
  # ID: 存储每个参与者分配到的 IMP 序号，初始为 NA
  ID <- matrix(NA, N, 1)
  
  for (i in 1:N) {
    # 修正海明距离计算：使用 sweep 安全地计算 x[i, ] 与所有 IRP 的距离
    # diff_matrix <- sweep(IRP, MARGIN = 2, STATS = x[i, ], FUN = "-")
    # distances[i, ] <- rowSums(abs(diff_matrix))
    distances[i, ] <- apply(IRP, 1, function(irp_row) {
      sum(abs(x[i, ] - irp_row))
    })
    # 找出最小距离的位置
    min_dist <- min(distances[i, ])
    min_posi <- which(distances[i, ] == min_dist)
    
    # 4. 修改后的逻辑：处理多重最小值 (Ties)
    if (length(min_posi) > 1) {
      # 存在多重最小值时，直接跳过并保留 NA (因为 ID 初始化为 NA，此处无需额外操作)
      ID[i] <- NA # 显式赋值 NA，以防万一
    } else {
      # 只有一个最小距离，直接分配 IMP 序号
      ID[i] <- min_posi
    }
  }
  # 返回对应的 IMP (掌握模式) 矩阵，对于 NA 的 ID，IMP[NA, ] 将返回 NA 组成的行
  # 为了避免 IMP[NA, ] 带来的复杂性，我们返回一个 N x K 的结果矩阵
  # ID_valid <- ID[!is.na(ID)]
  # # 计算每类对应的行号
  # ID_groups <- lapply(1:NI, function(k) which(ID == k))
  ID_groups <- lapply(seq_len(NI), function(k) which(ID == k))
  # result_IMP <- matrix(NA, N, ncol(IMP))
  # 
  # # 仅对成功分配 ID 的参与者进行赋值
  # assigned_indices <- which(!is.na(ID))
  # if (length(assigned_indices) > 0) {
  #   # 确保 ID 是数值型才能正确索引
  #   result_IMP[assigned_indices, ] <- IMP[as.numeric(ID[assigned_indices]), ]
  # }
  return(ID_groups)
}


# 比率返回 不同Q 不同的海明距离
distinguish_Q <- function(LCA_result_people, HDD_factor) {
      lv_max <- c()
      for (j in 1:length(LCA_result_people)) {
        LCA_people_row <- unlist(LCA_result_people[j])
        pairing_all_lv <- sapply(HDD_factor,function(x) {
          length(intersect(LCA_people_row,x))/length(LCA_people_row)})
        lv <- max(pairing_all_lv)
        lv_max <- c(lv_max,lv)
      }
  return(lv_max)
}


match_rows_simple <- function(response_Initial) {
  J <- ncol(response_Initial) #题目个数
  same_matrix <- matrix(data = NA,J,J)#预设空矩阵是题与题之间的联系
  row.names(same_matrix) <- paste0("Q",1:J)#矩阵命名
  colnames(same_matrix) <- paste0("Q",1:J)
  for (i in 1:J) {
    for (j in 1:J) {
      same_matrix[i,j] <- sum(response_Initial[,i] == response_Initial[,j])
    }
  }
  diag(same_matrix) <- 0  #对角赋值为0
  # 对于每一列取出最大数的行名，即对题目答题模式的相似性，找到每个题目最相似的对应分类题目
  max_positions <- sapply(1:ncol(same_matrix), function(i) which.max(same_matrix[, i]))
  #sapply(1:ncol(same_matrix), function(i) rownames(same_matrix)[which.max(same_matrix[, i])])
  max_positions <- data.frame(max_positions)
  item <- c(paste0("Q",1:J))
  max_positions_finish <- cbind(item,max_positions)
  max_positions_finish$max_positions <- sprintf("Q%d", max_positions_finish$max_positions)
  mat <- max_positions_finish
  # 上述操作将每道题目对应最相似的题目提出来，看能够提取出对于题目的分类
  
  matched_groups <- list()
  distance <- list()
  df <- as.data.frame(mat)
  processed <- logical(nrow(df)) # 作为循环逻辑判别
  
  # 使用for循环遍历所有行
  for (current_row_index in 1:nrow(df)) {
    if (processed[current_row_index]) {
      next  # 如果当前行已经处理过，则跳过
    }
    current_group <- as.character(df[current_row_index, ])
    processed[current_row_index] <- TRUE   #索引变成true
    
    repeat {
      new_overlap_found <- FALSE # 刚开始没有重叠行
      for (i in 1:nrow(df)) {
        if (!processed[i] && any(current_group %in% as.character(df[i, ]))) {
          # 如果有重叠，上面进行逻辑判定后，合并行
          current_group <- unique(c(current_group, as.character(df[i, ])))
          # 标记该行为已处理
          processed[i] <- TRUE
          new_overlap_found <- TRUE  # 找到新重叠
        }
      }
      
      if (!new_overlap_found) {
        break
      }
    }
    matched_groups <- append(matched_groups, list(current_group))
  }
  return(matched_groups)
}
























